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The effects of polymer additives on Rayleigh-Taylor (RT) instability of immiscible fluids 
is investigated using the Oldroyd-B viscoelastic model. Analytic results obtained ex- 
ploiting the phase-field approach show that in polymer solution the growth rate of the 
instability speeds up with elasticity (but remains slower than in the pure solvent case). 
Numerical simulations of the viscoelastic binary fluid model confirm this picture. 
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1. Introduction 

Mixing of species (e.g. contaminants, tracers and particles) and thermodynamica l 
quantities (e.g. temperature) are dramatically influenced by fluid flows (jPimotakis 112005 ). 
Controlling the rate of mixing in a flow is an objective of paramount importance in many 
fields of science and technologies with wide-ranging consequences in industrial applica- 
tions (jWarnatz. Maas fc Dibble Il200lh . 

The difficulties of the problem come from the intricate nature of the u nderlyin g fluid 
flow, which involves many active nonlinear ly coupled degrees of freedom (jFrisch Ill995h . 
and on the poor comprehension of the way through which the fluid is coupled to the 
transported quantities. The problem is even more difficult when the transported quan- 
tity reacts back to the flow field t hus affecting its dynamics. An instance is provided by 
the heat transport in convection (jSiggia Ill994l ). 



Mixin g emerges as a final stage of successive hydrodynamic instabilities (jDrazin fc Reid I 
19811) eventually leading to a fully developed turbulent stage. The possibility of control- 
ling such instability mechanisms thus allows one to have a direct control on the mixing 
process. In some cases the challenge is to enhance the mixing process by stimulating 
the turbulence transition, in yet other cases the goal is to su ppress deleterious ins tabili- 
ties and the ensuing turbulence. Inertial confinement fusion ( Cook fc Zhoull2002f ) is an 
example whose success relies on the control of the famous Rayleigh-Taylor (RT) insta- 
bility occurring when a heavy, denser, fluid is accelerated into a lighter one. For a fluid 
in a gravitation al field, such instability was first described Lord Rayleigh in the 1880s 
(I Rayleigh 18831) an d later generalized to all accelerated fluids by Sir Geoffrey Taylor in 
1950 (lTavlorlll950l) . 

Our attention here is focused on RT instability with the aim of enhancing the pertur- 
bation growth-rate in its early stage of evolution. The idea is to inject polymers into the 
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fluid and to study on both analytical and numerical ground how the stability of the re- 
sulting viscoelastic fluid is modified. Similar problems were already investigated in more 
specific context, including RT instability of viscoelastic fluids with suspended particles 
in porous medium with a magnetic field (jSharma fc Rajput 1992f) and RT linear sta- 
bility analysis of viscoelastic drops in high-speed airstream (jjoseph. Beavers fe Funada 



2002). We also mention that the viscoelasticity is known to affect also other kin d of insta- 



bilities, inchjdjng_Saffinaii^ (Wilson 1ll990t ICoussot P.l ll999l), Faraday 

(|Wagner. Muller fc Knorr lll999t iMuller fc Zimmermann I Il999h. the stability of 



waves 



Kolmogorov flow (iBoffetta et al. 20051 ). Tavlor-Couette flow dLarson. Shaqfeh fc Muller 



1990; Groisman fc Steinberg Ill996l ). and Rayleigh-Benard problem (jVest fc Arpaci II 19691 : 
Sokolov fc Tanner Ill972h 



The paper is organized as follows. In Sec. 2 the basic equations ruling the viscoelastic 
immiscible RT system are introduced together with the phase-field approach. In Sec. 3 
the linear analysis is presented and the analytical results shown and discussed in Sec. 4. 
The resulting scenario is corroborated in Sec. 5 by means of direct numerical simulations 
of the original field equations. 



2. Governing equations 

The system we consider is composed of two incompressible fluids (labeled by 1 and 2) 
having different densities, p\ and P2 > Pi, and different dynamical viscosities, p\ and p,i, 
with the denser fluid placed above the less dense one. For more generality, the two fluids 
are supposed to be immiscible so that the surface tension on the interface separating the 
two fluids will be explicitly taken into account. 

The effects of polymer additives is here studied w ithin the framework of the Oldroyd-B 



1 ne ertects or polym er additives is ncre studied wi 
model (|Oldrovdlll950t lHinchlll977t iBird et a/.l[i"987T) 



In this model polymers are treated 
as elastic dumbbells, i.e. identical pairs of microscopic beads connected by harmonic 
springs. Their concentration is supposed to be low enough to neglect polymer-polymer 
interactions. The polymer solution is then regarded as a continuous medium, in which 
the reaction of polymers on the f low is described as an elastic contribution to the total 



stress tensor of the fluid (see e.g. IBird et aLlll987t ). 



In order to describe the mixing process of the resulting viscoelastic immiscible flu- 
ids we fol l ow th e phase-field approach ( for a general description of the method see, 
e.g. iBravl (l2002h: ICahn fc Hilliardl (ll958Tl. a n d for application to multi p hase flows see , 
e.gjBadalassi. Ceniceros fc Baneriee I (j2003l) ; IPing. Spelt fc Shu I (|2007l) ; iMorrol ([20071 ); 
Celani et all (|2009h ). Here, we only recall that the basic idea of the method is to treat 
the interface between two immiscible fluids as a thin mixing layer across which physical 
properties vary steeply but continuously. The evolution of the mixing la yer is ruled by an 
order parameter (the phase field) that obeys a Cahn-Hilliard equation (|Cahn fc Hilliard 
19581) . One of the advantage of the method is that the boundary conditions at the fluids 
interface need not to be specified being encoded in the governing equations. From a nu- 
merical point of view, the method permits to avoid a direct tracking of the interface and 
easily produces the correct interfacial tension from the mixing-layer free energy. 

To be more specific, the evolution of the viscoelastic binary fluid is described by the 
system of differential equations 



p (d t v + v ■ dv) = -dp + d ■ (2/ie) + Ap a g<j) - cj>dM 



2pri 



8 • (<r - I) 



(2.1) 



d t (f> + v ■ d(f> = jd 2 M 



(2.2) 
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d t cr + v d(T = (3v) T cr + er dv -(cr - 1) . (2.3) 

r 

Eq. (|2.ip is the usual Boussinesq Navier-Stokes equation ( Kundu fc Cohen 200 ll ) with 
two additional stress contributions. The first one arises at the interface where the effect 



of surface tension enters into play (jBravl 120021 : 1 Yue et al. 1 12004; Be rti et al. 1120051) . the 
last term represents the polymer back-reaction to the flow field (jBird et al. 1987[ ). 
In (|2.ip . we have defined po — (pi + P2)/2, g is the gravitational acceleration pointing 
along the y-axis, A = (p2 — Pi)/{p2 + pi) is the Atwood number, = (diVj + djVi) /2 
is the rate of strain tensor an d p = p((f>) is the dynamical viscosity field parametrically 
defined as (|Liu fc Shen 1 120031) 



f 



2pi 



+ 



2p 2 



(2.4) 



4> being the phase field governed by (|2.2| . The phase field <ft is representative of density 
fluctuations and we take cf> = 1 in the regions of density pi and (f> = — 1 in those of 
density 02 > pi- cr = '^ffl is the polymer conformation tensor, _R being the end-to-end 
polymer vector (Rq is the polymer length at equilibrium) , the parameter r\ is proportional 
to polymer concentration and r = r((t>) is the (slowes t) polymer relaxation time which, 
according to the Zimm model ( Doi fc Edwards lll986l ). is assumed to be proportional to 
the viscosity p (therefore we have r = 7i for cj> = 1 and r = T2 for <fi = —I with /i(0)/r(0) 
constant). Finally, 7 is the mobility and M. is the chemical potential defined in terms 
of th e Ginzburg-Landau free energy T as (jCahn fc Hilliardl[f95i lBravll2002t lYue et al. 
120041) 

M = 6 -^ and ^[0] =A^daj Q|90| 2 + ^(<^)^ . (2.5) 

where Q is the region of space occupied by the system, A is the magnitude of the free- 
energy and the potential V(4>) is 



v{<t>) 



1 



I) 2 



(2.6) 



where e is the capillary width, representative of the interface thickness. 



by 



The unstable equilibrium state with heavy fluid placed on the top of light fluid is given 



v = , <p(y) = — tanh ( 



and 



(2.7) 



corresponding to a planar interface of width e with polymers having t heir equilibrium 
lengt h Rq. In this case, the surface tension, S, is given by (see, for example. lLandau fc Lifshitz 
2000): 

S = \f~dy Q|d</>| 2 + V(0))=^. (2.8) 

The sharp-interface limit is obtained b y taking the A and e to zero, keeping S fixed to 
the value prescribed by surface tension (|Liu fc Shen 1120031 ). 



3. Linear stability analysis 

Let us now suppose to impose a small perturbation on the interface separating the 
two fluids. Such perturbation will displace the phase field from the previous equilibrium 
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configuration, which minimizes the free energy ()2.5|) to a new configuration for which, in 
general, M 7^ 0. We want to determine how the perturbation evolves in time. 

Focusing on the two-dimensional case (corresponding to translational invariant per- 
turbations along the z direction), let us denote by h(x,t) the perturbation imposed to 
the planar interface y = in a way that we can rewrite the phase-field 4> as: 



/ 



y-h(x,t) 
eV2 



(3.1) 



where h can be larger than e, yet it has to be smaller than the scale of variation of h 
(small amplitudes). In this limit we assume the interface to be locally in equilibrium, 
i.e. d 2 f/dy 2 — V'(f), and thus f(y) — — tanh(y) and therefore M = — A|^j (' denotes 
derivative with respect to the argument). 



Linearizing the momentum equation for small interface velocity we have 
PodtUy = -d y p - 4>d y M - Agp <j) + ^P-dia i2 + p(d 2 + d 2 ) v y + 2{d y v y )d v f 
Integrating on the vertical direction and using derivations by parts one gets 



pod t q = S 



d 2 h 
dx 2 



2Agpoh 



2/i?7, 



where we have defined 
/■+00 

Q = 



p 



d 2 d 2 , , 



v y dy 



Q 



£ = 



(3.3) 



d x ai2 dy , (3.4) 



and we have used the relations / {f'fdy = 2 v / 2/(3e), / ff'dy = 0, / fdy = 2h. 

Note that, unlike what happens in the inviscid case, Eq. (|3.3[) does not involve solely 
the field q y but also second-order derivatives of v y . In order to close the equation, let us 
resort to a potential-flow description. The idea is to ev aluate Q for a potential flow v y 
and then to plug Q = Q pot into (|3.3p (jMikaelian II1993T ). The approximation is justified 
when viscosity is sufficiently small and its effects are confined in a narrow region around 
the interface. Because for a potential flow d 2 v = we have 

r0 



Q 



pot 



p- 



8 Uy 

dx 2 



dy 



/'- 



d m„ 



dy = (pi + p 2 



d q 



Substituting in (|3.3[) and defining v = (p\ 



dx 2 dx 2 
P2)/(2po) one finally obtains 



S d 2 h 
po ox 1 



2Aqh+ 2 ^ + 2^. 

Tpo OX* 



(3.5) 



(3.6) 



Let us now exploit the equation (|2.2|) for the phase field to relate q y to h. For small 
amplitudes, we have: 



and therefore, from 



d 2 M 



-f'd t h 



X 



, d 4 h 

dx 1 ' 2e 



l_ fll ,cPh 

2 Q x 2 



7A 

e 



f'dth- 



;f"d 2 x h 



Integrating over y, observing that l/(2v / 2e)/' approaches S(y — h) as e 
the limit of sharp interface (7 A — > 0) one obtains 

8th = v y (x, h(t, x), t) = vf 1 ^ (x,t) . 



(3.7) 

(3.8) 
and using 

(3.9) 
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The equation for the perturbation ai 2 of the conformation tensor is obtained by lin- 
earizing (|2.3|) around a a p = 5 a p 

dt(Ji2 — d x v y + dyV x - -CT12 (3.10) 

from which, exploiting incompressibility, we obtain 

2 1 
d t d x a 12 = {dl - dl)v y d x a 12 - 2a 12 d x - . (3.11) 

T T 

For small amplitude perturbations the last term, which is proportional to (Ji 2 d x 4>, can 
be neglected at the leading order. Integrating over y and using again the potential flow 
approximation one ends up with 

d t E = 2d 2 x q - ?E - (1 - 1) f dy<j)d x a 12 . (3.12) 

t n t 2 J 

where we have introduced f = 2T\T 2 j (n + t<£). 

In conclusion, we have the following set of equations (in the (x, t) variables) for the 
linear evolution of the Rayleigh-Taylor instability in a viscoelastic flow 



dth = Vy 



(int) 



d t q = f a dlh + 2Agh+^Y> + 2vdlq (3.13) 
where c = 4^i/x 2 /(jUi + H2) 2 < 1- 

4. Potential flow closure for the interface velocity 

The set of equations (|3.13|) is not closed because of the presence of the interface velocity 
v y mt ^ and of the integral term in the equation for E. In order to close the system we exploit 
again the potential flow approximation for which v y — d y ?p. 

Taking into account the boundary condition for y — > 00, the potential can be written 
(e.g. for y > 0) as 



ip(x,y,t)= e- ky+lkx i;{k,t)dk + c.c. (4.1) 
Jo 

where denotes the Fourier transform, and therefore 

/>oo 

v y {x,y,t) =- ke- ky+lkx *P(k,t)dk + c.c. (4.2) 



oo 

ikx _ 



q{x,t) = -2 e lkx *l)(k,t)dk + c.c. (4.3) 
and taking a flat interface, y = 0, at the leading order 

/>oo 

w (mt) (x, t) = - ke lkx ^P(k, t)dk + c.c. (4.4) 
Jo 

Assuming consistently that also 

/>oo 

a 12 (x,y,t)= / e- k y+ tkx a 12 (k,t)dk + c.c, (4.5) 
Jo 

in the limit of small amplitudes one has J dy4>d x a\ 2 = and the set of equation (|3.13[) 
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for the Fourier coefficients becomes 
d t h = \q 

3 t q = ~-j-k 2 h + 2Agh+^f-t-2vk 2 q (4.6) 
d t ± = -2k 2 q-\t. 

Restricting first to the case without polymers (?) = 0), the growth rate a at of the pertur- 
bation is obtained by looking for a solution of the form h ~ e aNt which gives 



a N = -vk 2 + ^Jlo 2 + (vk 2 ) 2 (4.7) 

where it has been defined 



w = J Agk - ^-P . (4.8) 

The expressio n (|4.8|) is the well-know n growth rate for a Newtonian fluid in the limit of 
zero viscosity (jChandrasekhar II1961J ). while (|4.7p is a known upper bound to the growth 



rate for the case with finite viscosity IjMenikoff et al. 1119771 ). 



Let us now consider the case with polymers, i.e. rj > 0. The growth rate a is given by 
the solution of 

(a?) 3 + 2(af) 2 (l + vk 2 f) + a? [4v(l + crf)k 2 f - uj 2 ? 2 ] - 2uj 2 ? 2 = . (4.9) 

The general solution is rather complicated and not very enlightening. In the limit of stiff 
polymers, f — * 0, one gets 

a = lim a = -v{\ + cq)k 2 + \] uj 2 + [v{l + crj)k 2 } 2 . (4-10) 

T— >0 



Comparing with (|4.7j) one sees that in this limit polymers simply renormalize solvent 
viscosity. This result is in agreement with the phenomcnological definitio n of c rj as the 



zero-shear polymer contribution to the total viscosity of the mixture (jVirk 1119751 ) . There- 
fore, in order to quantify the effects of elasticity on RT instability, the growth rate for 
viscoelastic cases at finite f has to be compared with the Newtonian case with renormal- 
ized viscosity v(l + erf). 

Another interesting limit is f — > oo. In this case from (|4.9|l one easily obtains that the 
growth rate coincides with that of the pure solvent (|4.7|) . i.e. = aN- The physical 
interpretation is that in the limit f — > oo and at finite time for which polymer elongation 
is finite, the last term in (|2.I[) vanishes and one recovers the Newtonian case without 
polymers (i.e. rj = 0). Of course, this does not mean that in general polymer effects for 
high elasticity disappear. Indeed in the long-time limit polymer elongation is able to 
compensate the 1/t coefficient and in the late, non-linear stages, one expects to observe 
strong polymer effects at high elasticity. 

From equation (|4.9|) one can easily show (using implicit differentiation) that a(f) is a 
monotonic function and, because OLoo ^ c^O; we have that instability rate grows with the 
elasticity, or the Deborah number, here defined as De = lot. 

The case of stable stratification, g — > —g, is obtained by to 2 — > — lo 2 neglecting surface 
tension. In this case (|4.9p has no solution for positive a, therefore polymers alone cannot 
induce instabilities in a stably stratified fluid. 



5. Numerical results 



The analytical results obtained in the previous Sections are not exact as they are based 
on a closure obtained from the potential flow approximation. While this approximation 
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is consistent for the inviscid limit v = (where it gives the correct result (|4.8|) for a 
Newtonian fluid) for finite viscosity we have sh own that it gives a kno wn upper bound 
to the actual growth rate of the perturbation ( Menikoff et al. 1977[ ) (this is because 
the potential flow approximation underestimates the role of viscosity which reduces the 
instability). Nonetheless, in the case of Newtonian fluid this upper bound is known to 
be a good a pproximation of the a ctual value of the growth rate measured in numerical 



JOT 

simulations (|Menikoff et al. Ill977| ). Because both r — ► and r — * oo limits correspond to 



Newtonian fluids, we expect that also in the viscoelastic case the potential flow description 
is a good approximation. 

To investigate this important point, we have performed a set of numerical simulations 
of the full model ()2.H2.3|) in the limit of constant viscosity and relaxation time (i.e. 
A*i = A*2, c = 1 and ri = T2 = f) in two dimensions by means of a standard, fully 
dealiased, pseudospectral method on a square doubly periodic domain. The resolution of 
the simulations is 1024 x 1024 collocation points (a comparative run at double resolution 
did not show substantial modificat ions on the results). More details on t he numerical 
simula tion method can be found in lCelani. Mazzino fc Vozella (2006) and lCelani et al\ 
(120091) . 

The basic state corresponds to a zero velocity field, a hyperbolic-tangent profile for 
the phase field and an uniform distribution of polymers in equilibrium, according to 
(|2.7p . The interface of the basic state is perturbed with a sinusoidal wave at wavenumbcr 
k (corresponding to maximal instability for the linear analysis) of amplitude ho much 
smaller than the wavelength (kho = 0.05). 

The growth rate a of the perturbation is measured directly by fitting the height of 
the perturbed interface at different times with an exponential law. For given values of 
A g 7 S/po, v and rj, this procedure is repeated for different values of f at the maximal 
instability wavenumber k (which, for the range of parameters considered here, is always 
k = 1, i.e. it is not affected by elasticity). Figure Q] shows the results for two sets of runs 
at different values of rj and v. As discussed above, we find that the theoretical prediction 
given by (|4.9| is indeed an upper bound for the actual growth rate of the perturbation. 
Nevertheless, the bound gives grow rates which are quite close to the numerical estimated 
values (the error is of the order of 10%). The error is smaller for t he runs having a larger 
value of rj and u, as was already discussed bv lCelani et al. I (l2009h . 

Both theoretical and numerical results show that the effect of polymers is to increase 
the perturbation growth-rate, a grows with the elasticity and saturates for sufficiently 
large value of De. 



6. Conclusions and perspectives 

We investigated the role of polymers on the linear phase of the Rayleigh-Taylor in- 
stability in an Oldroyd-B viscoelastic model. In the limit of vanishing Deborah number 
(i.e. vanishing polymer relaxation time) we recover a known upper bound for the growth 
rate of the perturbation in a viscous Newtonian fluid with modified viscosity. For finite 
elasticity, the growth rate is found to increase monotonically with the Deborah number 
reaching the solvent limit for high Deborah numbers. Our findings are corroborated by 
a set of direct numerical simulations on the viscoelastic Boussinesq Oldroyd-B model. 

Our analysis has been confined to the linear phase of the perturbation evolution. 
When the perturbation amplitude becomes sufficiently large, nonlinear effects enter 
into play and a fully developed turbulent regime rap idly sets in (|Cabot fc Cook 1 120061: 
Vladimirova fc Chertkov I [20091 : iBoffetta et al. Il2~009f ). In the turbulent stage we expect 



more dramatic effects of polymers. In turbulent flows, a spectacular consequence of vis- 
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Figure 1. The perturbation growth-rate a normalized with the inviscid growth rate lo (|4.8[) as 
a function of the Deborah number De = lot. Points are the results of numerical simulations of 
the full set of equations (|2.H2.3|I . lines represent the theoretical predictions obtained from (|4.9[) . 
The values of parameters are: c = 1, k = 1, Ag — 0.31, S/po — 0.019 and r\ = 0.3, v = 0.3 
(upper points and line) and n = 0.5, v = 0.6 (lower points and line). 



coelasticity induced by polymers is the drag reduction effect: addition of minute amounts 
(a few tenths of p. p.m. in weight) of long-chain soluble polymers to water leads to a strong 
reduction (u p to 80%) of the power n ecessary to maintain a given throughput in a chan- 
nel (see e.g. iToms" 19491 Virk I [l975l) . We conjecture that a similar phenomenon might 
arise also in the present context. Heuristically, the RT system can indeed be assimilated 
to a channel inside which vertical motion of thermal plumes is maintained by the avail- 
able potential energy. This analogy suggests the possibility to observe in the viscoelastic 
RT system a "drag" reduction (or mixing enhancement) phenomenon, i.e. an increase of 
the velocity of thermal plumes with respect to the Newtonian case. Whether or not this 
picture does apply to the fully developed turbulence regime is left for future research. 
We thank anonymous Referees for useful remarks. 
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